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The slow dynamics of nearly stationary patterns in a FitzHugh-Nagumo model are studied using a phase dy- 
namics approach. A Cross-Newell phase equation describing slow and weak modulations of periodic stationary 
solutions is derived. The derivation applies to the bistable, excitable, and the Turing unstable regimes. In the 
bistable case stability thresholds are obtained for the Eckhaus and the zigzag instabilities and for the transition 
to traveling waves. Neutral stability curves demonstrate the destabilization of stationary planar patterns at low 
wavenumbers to zigzag and traveling modes. Numerical solutions of the model system support the theoretical 
findings. 




I. INTRODUCTION 

Studies of stationary patterns in activator-inhibitor systems 
have focused primarily on localized structures such as pulses 
and spots in excitable and bistable media |[ P 
and periodic patterns near a Turing bifurcation 
Localized structures have instabilities to traveling patterns, 
breathing motion, and transverse deformations |0, |l2|, |l3| [ij] . 
Periodic patterns have been analyzed near the onset of a Tur- 
ing instability and also near the codimension-two point of a 
Turing instability and a Hopf bifurcation [|l5[ |l^ [l|, [l|] . 
But very few studies have explored instabilities of periodic 
(nonlocal) stationary patterns in excitable and bistable media, 
or of periodic stationary patterns far beyond the Turing in- 
stability [ pT| , p2| . The latter case includes pattern formation 
studies on the CIMA chemical reaction [[T7|,]23|, p4{ |. 

In this paper we study instabilities of stationary periodic 
patterns by deriving a Cross-Newell phase equation [ pO[ p5| , 
Pq ]. The derivation is not restricted to the immediate neigh- 
borhood of a Turing instability and applies to periodic patterns 
with space-scale separation that arise far from onset or in ex- 
citable and bistable media. The Cross-Newell equation was 
originally derived in the context of fluid dynamics and has re- 
cently been applied in a laser system p7[]. 

We choose to study the FitzHugh-Nagumo (FHN) equa- 
tions, a canonical model for activator-inhibitor systems. 
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'dt 
dv 
dt 



u — — V + V^u , 

e{u — aiv — flo) + (5V^w . 
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an excitable medium, a bistable medium, or a system with 
a Turing instability [pj|]. All three cases support stationary 
periodic solutions for S sufficiently large. 

In Section ^ we derive a phase equation describing weak 
modulations of periodic stripe pattern in the FHN model. In 
Section ^ we evaluate stability thresholds for Eckhaus and 
zig-zag instabilities and for a transition from stationary to 
traveling patterns. These thresholds suggest a number of spa- 
tial or spatio-temporal behaviors which we test in Section |^ 
with numerical solutions of Eqs. (IT]). 



II. THE PHASE EQUATION 

Let uo{9; k) = UQ{e + 27r; k), vo{9; k) vo{9 + 27r; k) 
be a stationary periodic solution of Eqs. (|l]) with phase 9 and 
wavenumber k. We consider weak spatial modulations of this 
periodic pattern and assume that those modulations have a 
length scale L that is much larger than the wavelength 1/k. 
The ratio of the length scales A = 1/ (kL) can then be used as 
a small parameter to write modulated solutions as an asymp- 
totic expansion about the periodic solution 



v( 



, R, T) = ua{9; k) + \ui{9, R, T) + \^U2{9, R, T) + ... 
I R, T) = vt){9; k) + \vi{9, R, T) + \^V2{9, R, T) + ...(2) 



Here, u is the activator and v the inhibitor. The parameters oq 
and ai can be chosen so that the FHN model dlj) represents 



where R = Ar and T = X^t are slow space and time vari- 
ables. The phase 9 in Eq. (^ is an undetermined function of 
space and time and fc = |k| = |V0| is the local wavenumber 
Our objective is to derive an equation for the slow phase 

e(R,r) := A0(r,R,r). 
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In terms of this phase the local wavevector is 

k(R,T) - v^e. 

Inserting the expansions (^) in Eqs. ([T]) we find at order unity 



uo-Uq-vo + 



e{uo - aivo - ao) + ok 



= 0, 
= 0, 



where /c" = k • k. At order A 
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Projecting the right hand side of (Q) onto {dgup, —e^^8gVo), 
the solution of the adjoint problem, produces the phase equa- 
tion 

r||=-VH(ki?), 

where 

T = < {8gUof > -e-i < {devo)^ > , (6) 
B = - < {8eUof > +Se-^ < {dgvof > ■ (7) 

In these equations < (.) >:= ^ J^^{.)d9. 

The quantities B and r contain information about vari- 
ous instabilities of the periodic stripe pattern. The condition 
■^[kB{k)] — implies the onset of an Eckhaus instability 



and the condition B — the onset of a zigzag instability | 
In Appendix A we show that the condition r = indicates a 
transition to traveling waves. 

To implement these conditions we need to solve Eqs. (||) for 
the periodic solution {uq, vq). For parameter values that sat- 
isfy e/S := fi ^ 1 an approximate solution can be computed 
as shown in Appendix B. Using this solution to calculate r and 
B, as shown in Appendix C, gives the following expressions: 



B = 



2V2 

3nk qirkri 
2^/2 v_ 



3nk qirky^ 



A-+A+ 



/3(A_)7(A_,A+), 
/3(A_)7(A_,A+), 



27r^ 



(8) 



(9) 



(10) 



where ^ = e/5,r] = Ve5, v± = (±1 — ao)/^^, = ai + 1/2, 
/3(a;) = cothqx — cschqa; , (11) 

and 

7(A_,A+) = -1 + i(l + ao)gA_cschgA„ 

+ ^(l-ao)'7A+cschgA+. (12) 

The quantities A+ and A_ denote the widths of domains with 
high and low values of u and v, respectively. The width is 
measured with respect to the spatial coordinate z = ^9 (see 
Appendix B). Given k, Eqs. and (|o|) can be solved for 
A+(fc) and A_(fc). Using these solutions in Eq. (^ graphs of 
r and kB as functions of k can be produced. 



III. STABILITY THRESHOLDS 



Explicit forms for r(fc) and B{k) are available in the sym- 

fc 



metric case, ao — 0, where A+ = A_ = — ir-: 



r(fc) 
B{k) = 



1 



where 77c 



Trkrjcq^ 
1 

nkricq^ 
and 



1 - ^f{7rq^/k) 
V 

1 - ^f{7Tq^/k) 



f{x) ~ (1 — X cscha;) (coth 2; — cscha;) . 



(13) 



(14) 



Figure |l] shows graphs of r(fc) and kB{k) for a bistable 
medium obtained with Eqs. (|3]) (thick lines) and with Eqs. (||) 
and (^ using numerically calculated solutions uq, vq (circles). 
A very good agreement is obtained within the validity range 
of the analysis, k 0{^) < 1. For fc the de- 

viations become large. In particular the minimum of kB{k) 
which designates the Eckhaus instability threshold, is not re- 
produced by the analytical form (pj|). 

The instability to traveling waves occurs at t = or at 



e = vif\T:q^/k)5-' . 
The zigzag instability occurs at B = or at 

e = r,lf{^q^/k)5. 



(15) 



(16) 



The condition -^{kB) 
comes 



for the Eckhaus instability be- 



^1 



0. 



Consider first the limit fc ^ in which the periodic pattern 
approaches an array of isolated front structures. In this limit 
fi'Kq-JJl/k) 1 and the condition for the onset of traveling 
waves becomes e — rjlS^^. This is precisely the nonequilib- 
rium Ising-Bloch (NIB) bifurcation point, where a stationary 
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FIG. 1: Typical functions T{k) and kB{k) for a bistable medium. 
The curves represent the functions of Eqs. Jl3). The circles are 
numerically computed solutions using Eqs. (^and The point 
kB — indicates the boundary between stable stationary stripes and 
zigzag patterns. At r = the pattern becomes unstable to traveling 
waves. Parameters: ai = 4, ao — 0, e — 0.001, 5 — 2.0. 



front loses stability to a pair of counter-propagating fronts. 
The condition for the zigzag instability becomes e — rj^ S. 
This is the threshold for the transverse front instability [p8|]. 

The neutral stability curves for a bistable medium corre- 
sponding to Eqs. (|l5]) and ( [1^ ) are shown in Figs. ^ and ^ 
for fixed S and e, respectively. They imply that high wavenum- 
ber stationary planar patterns are stabilized against zigzag and 
traveling wave instabilities. Notice that for 5=1 the neutral 
stability curves r = and B = coincide (see Eqs. ( |l5[ ) and 
(16) or Fig. ^). For S > I, upon decreasing the wavenum- 
ber at constant e, a high wavenumber pattern is destabilized to 
a zigzag pattern, whereas for 6 < 1 the destabilization is to 
traveling waves. Similar neutral stability curves are found for 
the nonsymmetric case, ao ^ 0, for excitable media and for 
systems (far) beyond the Turing instability. 



IV. COMPARISONS WITH NUMERICAL SOLUTIONS 

We have computed numerical solutions of Eqs. (|l]) to test 
the stabilization of zigzag and traveling-wave instabilities at 
high wavenumbers. Figure |] shows a low wavenumber zigzag 
pattern and a high wavenumber planar pattern computed for 
the same parameter values. This behavior is well known in 
other contexts [^^. The zigzag instability is a mechanism b_ 
which the system locally increases the wavenumber Figure 
shows coexistence of a low wavenumber traveling wave and 
a high wavenumber stationary pattern. These numerical re- 
sults are for a bistable system but similar results are found for 
excitable and Turing unstable systems. Coexistence of sta- 
tionary and traveling waves has been found in experiments on 
the CIMA reaction QlTi e3]. and analyzed using different the- 
oretical approaches pi|, |T]]. 

We have also tested the condition for the Eckhaus insta- 
bility in a bistable system using numerical computations of t 




FIG. 2: The neutral stability boundaries for the zigzag instability 
(B = 0; thick, solid curve) and traveling wave instability (r = 0; 
thin, dashed curve). To the left of the B — curve planar periodic 
patterns are unstable to zigzag patterns. To the left of the t — curve 
planar periodic patterns are unstable to traveling waves. Parameters: 
ai = 2, ao = 0, e = 0.01, 5 = 2. 



and B. Choosing wavenumbers k > kc where kc corresponds 
to the minimum of kB, we found initial periodic patterns ei- 
ther collapse to uniform states or to a lower wavenumber pat- 
tern through phase slips. Fig. ^ demonstrates these two cases. 
Similar conclusions hold for excitable systems. An unstable 
Turing pattern, on the other hand, always converges to a lower 
wavenumber pattern since the single uniform state is unstable. 



V. CONCLUSION 

We have shown that the Cross-Newell phase equation pro- 
vides a powerful tool for studying instabilities of stationary 
periodic patterns in activator-inhibitor systems. The equation 
contains information not only on the Eckhaus and zigzag in- 
stabilities, but also on the destabilization of stationary peri- 
odic patterns to traveling waves. The same equation applies 
to bistable, excitable, and Turing unstable systems. Equations 
of that kind should prove useful in identifying parameters and 
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FIG. 3; Coexistence of zigzag and planar patterns. The dark ar- 
eas indicate regions of it > and the light regions u < 0. At 
higher wavenumbers (top) the planar stripe solution is stable. At low 
wavenumbers (bottom) the planar solution is unstable and forms a 
zigzag pattern. Parameters: al = 2, ao = 0, e — 0.05, 5 = 2. 



initial conditions where zigzag and Eckhaus instabilities cou- 
ple to traveling wave modes. Such coupling may lead to com- 
plex spatiotemporal behavior analogous to the coupling of the 
NIB front bifurcation to a transverse front instability | p^ , p2| ] 
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FIG. 4: Coexistence of traveling waves and stationary waves. 
At high wavenumber (top) the patterns are stationary and at low 
wavenumber patterns (bottom) they travel. Parameters: al = 2, 
ao = 0, e = 0.03, (5 = 2. 



APPENDIX A: THE MEANING OF r = 

We show here that the condition t = defines the critical 
value of e at which traveling solutions bifurcate from the sta- 
tionary solution. We look for traveling solutions u{9),v{9) of 
Eqs. (jl]), where 9 = kx — cut, that bifurcate from the stationary 
solution branch a; = at some e = Cc- Near the bifurcation 
where cj <C 1 we can expand the traveling solutions as power 
series in uj around the stationary solution uo,vo: 

u{9) = uo{9) + ujui{9) + ... , 
v{9) = vo{9)+ujvi{9) + ..., (Al) 
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Projecting the right hand side onto (uq, — ^v'q) gives 
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FIG. 5: Time evolution of a periodic pattern in the region of Eck- 
iiaus instability. The high wavenumber pattern is unstable and either 
converges to one of the uniform states (top) or a lower wavenumber 
pattern (bottom). Parameters: al = 2, ao — 0, S — 2; top: e = 0.01, 
bottom: e — 0.1. 



Expanding e as 



e = ec + f-iLO - 



(A2) 



and using these expansions in Eqs. ^ we find at order 



UJ 



CcUi + (sk'^-^ - eai ) = -^o - ei(uo - aivo - ao) 



<v'o'> 



(A3) 



where we used Eq. ( pb| ) and switched to the notation of a 
prime for the derivative with respe ct to the single argument 
9. Using the definition (ph of r and ( A3 ) we find 



(A4) 



Thus, T — implies e = ec or the onset of traveling solutions. 

APPENDIX B: APPROXIMATE STATIONARY SOLUTION 

For /i = e/(5 <C 1 a singular perturbation approach can 
be used to approximate the stationary solution uq{9),vo{0). 
Rescahng the space coordinate as 2; = -^9, Eqs. (||) become 

uo - Uo - vq + fm'(^ = , 
uq - aiVQ - ao + = , 

where the prime denotes now the derivative with respect to z. 
Since the small parameter /i multiplies the second derivative 
term Uq, two types of spatial regions can be distinguished. 
Outer regions where uq{z) varies on a scale of order unity 
and the term /iUp is negligible, and inner regions where uo{z) 
varies on a very short scale of order y/JI and the term ^Uq 
cannot be neglected. In these regions, however, vq hardly 
changes. 

The analysis of the inner regions leads to the solutions 



Uq = ± tanh 



V2k' 



vo = 0. 



(Bl) 



These solution represent front structures separating two types 
of outer regions: domains of high activator values, u — 
u+{vq) ("up state"), and domains of low activator values 
u — u_(wo) ("down state"), where u±{vo) are the extreme 
0. We look for periodic stationary 



roots of Uq — Uq — Vq 

solutions with wavelength A = A_ + A+, where A+ and A_ 
are the widths of up and down states respectively. Consider 
now a down state spanning the spatial range ~A_ < z < 
followed by an up state spanning the range < z < A+. The 
equations for v at these outer regions are 



g^(wo-w-) = 0, -A_<z<0, 



(B2) 



with the boundary conditions ?7o(— A_) = vq{0) ~ 0, and 

<-<z'(«o-t'+) =0, 0<z<A+, (B3) 

with the boundary conditions Wo (0) = uo(A+) = 0. In obtain- 
ing these equations we approximated u±{vo) — ±1 — uo/2. 
This approximation is particularly good for bistable media 
with ao small and ai relatively large. These values restrict 
Vq to a small range around vq = 0. For excitable media and 
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systems undergoing Turing instability, a large value of 5 might 
be needed to keep vq small. 

The solutions to Eqs. (B2) and (B3) are 



vo 



sinh qA 
for -A- < z < 0, and 

v+ 



sinh — sinh (/(z + A_)] + i;_ , (B4) 



■^0 



sinh qA^ 



sinhg(z — A_|_) — sinhgz] + -y+ , (B5) 



for < z < A+. To determine A± for a given A we match 
the derivatives of vq at the front positions 

v',{0+)=v'„{0-), v',{A+)^v'„{^A^). 

This leads to the relation 

v+P{qA+)+v^P{qA^) = 0, 

where f3{x) is given by Eq. (pi]). 

APPENDIX C: CALCULATION OF r AND B 

The quantities r and B are given by Eqs. and (Q). Con- 
sider first the integral 



1 

2^ 



It has a contribution from two inner regions at z = and 
z = A-)- where uq is given by Eq. (Bl), and a contribution 
from two outer regions, — A_ < z < and < z < A+ 
where uq = —l — vo/2 and uq = 1 — vo/2 with vq given by 



Eq. (U) and Eq. (§1), respectively. (Recall that z = ^0). 
The contribution from the two inner regions is 



27rfc2 



sech 



V2k 



de. 



V^TTk J- 



sech xdx 



2V2 



We have used here the fact that fc 0{y/JI) 1. The integral 
over a narrow inner region is transformed into an integral over 
a wide region after stretching the 6 variable to the x — 
variable. The contribution from the two outer regions is 



< u'q{9)'^ >oute 



8nk 



where we used in the two outer regions uo(z)' — 
Altogether, 

2V2 ^ 



37rfc 87rfc 




VQ{z)''^dz+ I vo{zy'^dz 



ivoizy. 



(CD 



A4 



Vo(z)'^dz+ / i;o(z)'^rfz 

-A_ Jo 

The second term on the right hand side of Eq. ( |Cl| ) is small 
(since ^/JI <^ 1) and will not contribute to the leading order 
forms of T and B. 

Consider now the integral 



The contribution from the inner regions to this integral is neg- 
ligible for /I <C 1. Thus 





vo{zy'^dz+ / voizy'^ dz\C2) 

-A_ Jo 



Using the solutions (B4) and (B5) in the integrals ( |C1[ ) 
and (C2) and using the expressions for r and B we obtain 
the expressions (p|). 
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